On the stationarity of linearly forced turbulence in finite domains 



O 



E. GravanifQ and E. AkylasQ 

Department of Civil Engineering and Geomatics, Cyprus University of Technology, P.O. Box 50329, 3603, Limassol, Cyprus 

(Dated: January 20, 2013) 

A simple scheme of forcing turbulence away from decay was introduced by Lundgren some time 
ago, the 'linear forcing', which amounts to a force term linear in the velocity field with a constant 
coefficient. The evolution of linearly forced turbulence towards a stationary final state, as indicated 
by direct numerical simulations (DNS), is examined from a theoretical point of view based on 
symmetry arguments. In order to follow closely the DNS the flow is assumed to live in a cubic 
domain with periodic boundary conditions. The simplicity of the linear forcing scheme allows 
one to re-write the problem as one of decaying turbulence with a decreasing viscosity. Scaling 
symmetry considerations suggest that the system evolves to a stationary state, evolution that may be 
understood as the gradual breaking of a larger approximate symmetry to a smaller exact symmetry. 
The same arguments show that the finiteness of the domain is intimately related to the evolution 
of the system to a stationary state at late times, as well as the consistency of this state with a high 
degree of isotropy imposed by the symmetries of the domain itself. The fluctuations observed in the 
DNS for all quantities in the stationary state can be associated with deviations from isotropy. Indeed, 
self-preserving isotropic turbulence models are used to study evolution from a direct dynamical point 
of view, emphasizing the naturalness of the Taylor microscale as a self-similarity scale in this system. 
In this context the stationary state emerges as a stable fixed point. Self-preservation seems to be 
the reason behind a noted similarity of the third order structure function between the linearly forced 
and freely decaying turbulence, where again the finiteness of the domain plays an significant role. 
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Maintaining a turbulent now in a more or less stationary 
state, for better statistics in experiment or convenience in 
O-itheoretical considerations, requires forcing the flow, that is 
feeding it energy which balances dissipation happening at 
the smallest scales. In numerical simulations of incompress- 
>- ible isotropic turbulent flows one usually solves the Navier- 
C*~) Stokes equations in a cubic box (with periodic boundary 
^1 conditions). For an account of Direct Numerical Simula- 
^ tion (DNS) methods see [J; for a recent review on the 
. current isotropic turbulence statistics from DNS see @ . In 
most cases forcing takes the form of a force term in wave 
number space (spectral space) which vanishes for all but 
the smaller wave numbers i.e. one feeds energy at the 
. . largest scales of the turbulent flow in the box. The general 
^ concept is that the details of the larger scales are model 
k>( dependent but the details of all other scales, that is those 
j_J yyhere some universal laws may hold, depend only on the 
intrinsic dynamics of the Navier-Stokes at least for high 
Reynolds numbers. Presumably, by forcing turbulence one 
achieves satisfactory results for given a resolution for higher 
Reynolds numbers than in the freely decaying turbulence. 

There have been developed various kinds of forcing 
schemes. The simpler ones fiddle in a suitable manner the 
magnitude of velocity field, or the total energy of the lower 
wave number modes, imitating an energy input in the larger 
scales These models can be regarded as essentially 

deterministic in the sense that that there is no additional 
randomness introduced in the problem. There are also de- 



terministic models which explicitly introduce a force term 
in the Navier-Stokes, whose details are eith er p ostulated or 
derived by a postulated auxiliary model |7HlC| . In stochas- 
tic forcing models [IlT[l3j the details of the force term are 
determined by additional random variables following pre- 
scribed stochastic processes. Each of those models suffers 
from one or more from a set of problems such as, excessive 
fluctuations around stationarity, relatively long relaxation 
period to stationarity, persistent anisotropy, excessive dis- 
tortion of large scale motions, introduction of irrelevant 
features in the description of turbulence. A useful compar- 
ative discussion between certain deterministic and stochas- 
tic models can be found i n (l4| . 

Lundgren proposed in [15| that we may simplify the de- 
terministic models to the bare minimum, in some sense, 
assuming that the usually velocity dependent force term is 
merely proportional to the velocity field for all positions 
x, or all wave numbers k, and all times: f = Au, where 
A is plainly a constant. The 'linear forcing' scheme was 
further studied in [l6[ and [Vq . Its simple force term Au 
has the same form in both the spectral and physical space. 
Thus, unlike other forcing schemes, it may be used equally 
easy in cases that need to be solved directly in the physical 
space with boundary conditions different than periodic 16] . 
That feature could prove useful. Additionally, although in 
the linear forcing the injection of energy into the flow is 
not restricted to the larger scales, this scheme performs de- 
cently, and in fact possibly better, in the region between 
the inertial ran ge a nd the integral scale than other forc- 
ing schemes in 15|. From the theoretical point of view 
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what matters most is that, unlike limited spectral band- 
width forcing schemes, linear forcing does not introduce an 
additional length scale in the problem at the level of the 
Navier-Stokes equations (a length scale outside the equa- 
tions is of course introduced by the boundary conditions). 
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FIG. 1: A typical evolution of the energy production rate (solid 
line), dissipation rate (dashed line) is shown in the figure (a) and 
Taylor microscale Reynolds number in figure (b). The parame- 
ters chosen are A = 1, box size / = 2ir and viscosity v — 0.1. 

The performance of the linear forcing scheme with re- 
spect to its convergence properties was studied in consid- 
erable detail in [16| and useful remarks have been made 
in [l4[. The clear conclusion is that linear forcing results 
in relatively large fluctuations in the stationary phase. In- 
deed, a typical evolution of the energy production rate 2AK 
(where K is the total kinetic energy per unit mass) , the dis- 
sipation rate e and the Taylor microscale Reynolds number 
Rexis shown in fig.Q] [The details of the DNS can be found 
in [l7j]. From the practical point of view this is a disad- 
vantage as it requires longer simulations in order to obtain 
good statistics. Also, the stationary state is reached af- 
ter a relatively long transient period [f| requiring even 
more computational time. On the other hand, linear forc- 
ing leads to quite controllable situations in the stationary 
state: Given the scales of the problem i.e. the rate A, the 
cubic box size I and the viscosity v, the facts of the station- 
ary state are predictable. The balance between the energy 
production and dissipation, 2AK — e, is indeed observed 
on the (time-) average validating the very concept of a sta- 
tionary state; the dissipation length L £ = (2K) 2 Je turns 
out to be equal to the box size I within few percent error 
in all cases [lg; the Reynolds number Re^ = K 2 /(ev) may 
be re- written as \AL\jv at the stationary state, should 
then be roughly equal to \ of the natural order of Rez, 



served [13] • For example, the Taylor microscale Reynolds 
number ReA = (^Re^)^ is expected to be roughly equal 
to 25.7 for the run shown in fig. [T] Indeed the average of 
the Re^ time-series in fig. [T] differs only by few percent from 
that estimate. 

Even if we take stationarity for granted, its character- 
istics i.e., the relatively large fluctuations and the 'pre- 
dictability' of quantities describing the state of turbulence, 
certainly call for understanding. On the other hand the 
very existence of a stationary state in this scheme is a fairly 
intriguing matter. The long-time effect of the energy pro- 
duction competing with dissipation is not a priori clear. 
From the dynamical point of view, it is clear that the dissi- 
pation term vXI 2 u becomes stronger than the force term Au 
at scales smaller than (v/A)* ~ Re^ 1 /, but it is not clear 
whether energy which is produced at all other scales up to 
I will be dissipated by an adequate rate at those smaller 
scales. 

We will approach the problem as follows. The relative 
simplicity of linear forcing allows us to study its late-time 
evolution employing scaling symmetry arguments to an ex- 
tent enjoyed possibly only in freely decaying turbulence; in 
fact as we shall show there is a relationship between lin- 
early forced and freely decaying turbulence. A quite par- 
allel discussion between them can be made. The sections 
HlUVl will be devoted in presenting these arguments. The 
predictability, as we called it above, of the stationary state, 
is enlightened through those symmetry arguments, essen- 
tially on the basis that there is no intrinsic large length 
scale in the dynamical equations apart from that intro- 
duced by the boundary conditions i.e. the finite size / of 
the domain. Then remains the question why the fluctua- 
tions observed in the stationary state, as seen e.g. in the 
fig. [U are so large. We shall argue, as analytically as we 
can, that the fluctuations can be associated with the devia- 
tions from isotropy accumulated by this forcing at all scales 
between the scale (v/A)? and the domain size I (unlike the 
limited bandwidth forcing schemes which feed anisotropy 
only at the domain size scale where isotropy is already bro- 
ken). The method we shall use is to reduce the dynamical 
problem to a two-equation model. As a cross check of our 
previous conclusions, the stationary state re-emerges as a 
stable fixed point of the evolution, a byproduct of which is 
that fluctuations tend to be suppressed as long as turbu- 
lence is isotropic. This part of our discussion is presented 
mostly in sections IVTl and lVlIl It is interesting to note that, 
from various aspects, linearly forced turbulence seems to be 
a natural context for the direct application of various ideas 
that have been developed in the study of freely decaying 
turbulence, in fact one may dare say, an even more natural 
context. 

In terms of equations, a linearly forced incompressible 
flow with zero mean flow velocity is described by the 
Navier- Stokes equation 



du 

~dt 



+ (u- V)u 



f . 



-Vp + iA7 2 u + Au, 



(1) 



in this problem, Al /v, in all cases, as it is indeed ob- where the incompressibility condition reads V • u = 0; the 
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velocity field is solenoidal. A is a positive constant with 
dimensions of inverse time. The term ^4u is a curious 'anti- 
drag' force on fluid particles. As already mentioned we 
impose periodic boundary conditions: u(x, y, z) = u(x, y + 
I, z) — u(x, y, z+l) = u(x+l, y, z). That is, the flow evolves 
within a cubic domain with side equal to I obeying the given 
conditions on its boundary. (We will often refer to the 
cubic domain simply as the 'box'.) As we shall emphasize 
later on, the term 'cubic domain' is slightly misleading due 
to the periodic boundary conditions: the flow essentially 
evolves in a boundary-less space of finite size. There are 
no walls anywhere, this is why we describe the domain as 
finite instead of bounded. The problem we are interested 
in to determine the late-time state of the turbulent flow 
governed by these equations and conditions. 

The present work is organized as follows. In sections [TT| 
and IIIII a reformulation of the problem and an associated 
scaling symmetry is presented. In section IIVI the impli- 
cations of the scaling symmetry and of the symmetries of 
the domain for the late-time behavior of the ensemble av- 
erage correlators are discussed. In section [V] we restrict 
ourselves to isotropic turbulence, to argue in a more de- 
tailed manner for the stationary state as the final phase 
of the linearly forced turbulence, as described by the ex- 
act ensemble average correlators and taking into account 
the effects of the finiteness of the domain. In section I VII 
the expected behavior of the actual observables in DNS 
i.e. the box-averaged correlators, is discussed in relation to 
the properties of the ensemble average correlators estab- 
lished in the previous sections. In section fVIII we combine 
the powerful condition of isotropy with the (by now es- 
tablished) existence of fluctuations around stationarity: A 
complete self-preserving isotropic turbulence model is ob- 
tained and applied to study the fate of fluctuations at scales 
in the flow where isotropy holds. We close with a few ad- 
ditional remarks and discussing certain open issues of the 
problem in the final section. 



II. UNFORCED TURBULENCE WITH 
DECAYING VISCOSITY 

We shall proceed as follows. Mathematically, we may re- 
write the problem as an equation for a new field u' which 
we require to satisfy 

^ + (u'.V)u' = — Vy + i/W, (2) 
of p 

with u' being related to u by 

u' = Fu, (3) 

where F is a function of time F — F(t') to be determined. 
V • u = implies that V ■ u' = 0. 
Substituting (0) to © we get 

1 dF 1 du , 11^,1 ,_ 2 

^^ U+ F^ + (U ' V)U = -^P VP+ F^ U - (4) 



This equation becomes identical to ([TJ upon setting 
1 dF 

Fdt' = dt, — — = -A, v' = Fv. (5) 
F 2 dt' W 

The transformation of pressure, p' = F 2 p, follows from its 
Laplace equation constraint and cannot be regarded as an 
independent condition. 

For constant A we have that 

F = e~ At , and At' + 1 = e At , (6) 

where we fix one integration constant by setting t' = 
corresponding to t = 0. It is left understood due to the 
freedom of the second integration consttant that we may 
shift t arbitrarily. 
Viscosity v' reads 



in the unphysical time coordinate t' . The problem has be- 
come unforced turbulence with decaying viscosity. The way 
it decays, oc is crucial in what follows. 

Note that everything can be transformed back to the 
initial variables at all times, except t = oo, or t' = oo. This 
is a singular point of the transformation as F vanishes and 
the two forms of the problem are not equivalent. That 
would be a relevant subtlety only if we had to deal with 
the actual limit t — > oo. We will not need such a limit 
anywhere in our analysis. 

It also worth to note that a transformation t — » t' can 
generate only a term which is linear in u, unless the trans- 
formation depends itself on the velocity field. Therefore 
the possibility of such a generating transformation is inti- 
mately related to linear forcing. 



III. TIME-TRANSLATION INVARIANCE AND 
AN EXACT SCALING SYMMETRY AT LATE 
TIMES 

Note that equation (UJ) does not explicitly depend on 
time, forcing being a function of the velocity field alone. 
Therefore, depending on boundary conditions, this equa- 
tion may admit non-trivial solutions which are static or 
independent of time in a certain sense. As we are inter- 
ested in fully developed turbulence, the time-dependence 
in question will apply to statistically defined quantities. 

The Navier-Stokes equations we arrived at by transform- 
ing to time t' in the previous section explicitly reads 

« + (rf .vy— iv^+^w, <8) 

where p, v and A are constants and V ■ u' = 0. Apart 
from the fairly peculiar time-dependence of viscosity, that 
is, the energy dissipation rate decreases with time, the form 
of this equation is more familiar than that of equation ([TJ . 
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Now, the time-translation invariance of equation ([T]) 
translates to an exact scaling symmetry of equation ([§]). 
Even by inspection one may verify that the transformation 

t' -> e a t' , u' -> e- a u' , (9) 

for any constant a is an exact symmetry of the previous 
equation (necessarily, p' — > e~ 2a p') for times t' ^ A -1 . 
Then the integration constants in the relation between t 
and t' are irrelevant. 

The origin of the late-time scaling symmetry is clear: 
Shifting t means rescaling t' , at least for times t' A^ 1 . 
Shifting t is a symmetry of equation (TTJ , therefore rescaling 
t' must be a symmetry of equation ([5]), as it is indeed the 
case. It is important to remember that the symmetry ^ 
respects the periodic boundary conditions on the field u', 
therefore it is an exact symmetry of the problem. 

We may now forget equation (fTJ) for a little while and fo- 
cus on the unforced turbulence described by (|5J). Its scaling 
symmetry will allow us to draw certain conclusions about 
the late time behavior of the system. 



IV. SCALING SYMMETRIES, ASYMPTOTIC 
BEHAVIOR AND ISOTROPY 

In order to get a first idea why the symmetry can be 
useful that way, note that the product t'u' is invariant un- 
der the scaling @. Consider an arbitrarily chosen moment 
of time t' and the velocity field u' at that moment, and 
another moment t' = e a t' when velocity is u'. Invariance 
means: t'u' = t' u' Q . Equivalently we may write 

u' = it u . (10) 

Now in general a symmetry transformation moves us 
around the space of solutions. That is, all the previous 
relation means, is that if there is a solution with velocity 
Uq at time t' then there is another solution with veloc- 
ity field u' at time t' . i.e. in general u' and u' need not 
necessarily correspond to the same initial conditions. 

On the other hand, the symmetry holds for large times 
t 1 and t' Q . Even if it did not, that would be a convenient 
choice for the following reason. The initial time t' = is 
pushed into the remote past, and the behavior (|10[) might 
then be an exact asymptotic result for a large class of so- 
lutions, meaning irrespectively of their initial conditions. 
That implies that t'u' — t' u' is an actual constant at each 
point r in space depending only on the parameters of the 
equation and the boundary conditions. 

The constant in question is a vector. To be more specific, 
recalling that u' satisfies V • u' = 0, we need a solenoidal 
vector field in steady state which does not depend on initial 
conditions i.e. it is unique. Such a field must respect the 
symmetries of the boundary conditions, that is the symme- 
tries of the cube. There is no such thing: solenoidal vector 
fields have closed integral curves which can always be re- 
versed by reflections. We deduce then that ([TO)) , as long as 



it is non-trivial, will always depend to some extend on t' 
i.e. on initial conditions. Thus it is not of much use in this 
form. 

Our reasoning can be used more effectively if it is applied 
in statistically defined quantities, that is, correlators of the 
velocity field. As mentioned already in the Introduction, 
from here and up until section fVTl we shall work with corre- 
lators defined as averages over a statistical ensemble. The 
statistical ensemble averages are independent of the initial 
conditions by their very definition: they are averages over 
the space of solutions. Of course in a problem on turbulence 
they certainly are the quantities of interest. The statistical 
ensemble averages will be denoted by an overbar. 

Then (fTOf holds trivially for no mean flow: u' = 0. 
Next one considers gen eral correlators of the velocity field, 
u' ix (ri,t[)u' i2 (r 2 , t' 2 ) ■ ■ •, and their derivatives. Consider lo- 
cal correlators i.e. all times and positions coincide. These 
are tensor fields Tj ij2 ... (r, t'). Let such a tensor field with n 
velocity field insertions in the correlation. Then the sym- 
metry © tells us, similarly to (TTU1) . that 



Now *o n 3ojija— ( r '*o) mus t be a constant at each point r 
in space. If not, then this quantity does depend on t' i.e. 
on the initial conditions. That means: this quantity is 
not well defined as a ensemble average i.e. mathematically 
does not exist and it must be defined in an approximate 
manner which does not possess the expected properties, or 
only approximately. The reason why this may happen is 
that the system has not reached a stage where ensemble 
averages are meaningful, a priori some kind of equilibrium 
is required. 

Now, same as with t' Q u' Q , most of these constant tensor 
fields must be zero by being inconsistent with the sym- 
metries of the cube (especially reflections) and the incom- 
pressibility condition. Certainly everything with at least 
one solenoidal index must vanish. This leaves us with 
the scalars, tensors manufactured out of them and Kro- 
necker delta, and correlators such as diUkdjUk with no free 
solenoidal index. 

In order to see how these statements are realized by an 
example, consider the correlator t^ u' Qi u' Qk , which is con- 
stant in time. Being constant in time means that it must 
respect the symmetries of the cubic domain: it must not 
change under reflections of the domain around planes of 
symmetry and rotations around axis of symmetry. One 
should recall that our correlators are ensemble averages 
over the whole of phase space, thus symmetries cannot take 
us to an other constant late-time solution: there is no other 
solution, or we have convergence problems in the very defi- 
nition of our averages. It is then easy to see that i 2 u' 0i u' ak 

must be equal to 5 t k i 2 u' 0l u' 0l (no sum) = |% u'^u^, 
i.e. essentially a scalar. Moreover, by the incompressibility 
condition V • u' = we see that the scalar itself must be 
constant in space. 

One should note that the situation resembles very much 
that of isotropic i.e. also homogeneous turbulence. There 
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is anisotropy allowed by the problem but it is much less 
than what would call anisotropy in general. Thus we will 
proceed by assuming isotropy and analyze what that im- 
plies; then, as isotropy cannot hold at scales comparable to 
the cubic box size I, the effects of the boundary eventually 
play a key role. This is done in the next section. We close 
this section by defining a few important scalars for the de- 
scription of turbulence, their symmetry and transformation 
properties and their expected late time behavior according 
to our arguments. 

The r.m.s. value q of the velocity and the dissipation rate 
e are defined by q 2 — u • u and e = v djUidjUi. Also by 
K = \q 2 we shall denote the total kinetic energy per unit 
mass. Similar expressions hold for the primed quantities. 

Under the symmetry ([9]) the quantities K 1 and e' trans- 
forms as 

K' -> e- 2a K' , e' -> e~ 3a e' , (12) 

where one should bear in mind that e' involves v' defined 
in (0 . Following again the reasoning given in the previous 
paragraphs we conclude that for large times t' the kinetic 
energy and dissipation rate should obey 

, constant , constant 



In order to see what this result means back in the vari- 
ables of the system ([TJ, we use (|3j and ([7j to obtain the 
transformations of K and e: 

K' = (At'+1)- 2 K, e' = (At' + iy 3 e. (14) 

The result is then that the kinetic energy and dissipation 
rate in the linearly forced turbulence should at late times 
become 

K = constant , e = constant . (15) 

Presumably, the dissipation length scale L e and the 
Reynolds number Rei, defined by 

a 3 K 2 
L e = 2-, ReL = — , (16) 

and transforming by 

Rc'l = Re L , L' E =L E , (17) 

should also reach constant values. That is, turbulence 
should get to what we have already called as the stationary 
state or phase. 

*** 

The arguments given above can be rephrased in the ac- 
tual time t and the variables of equation ([TJ as follows. We 
have already mentioned that shifting time t is a symmetry 
of equation ((TJ. That is, if u(t) is a solution of this equation 
then so is u(t+ At) for an arbitrary interval At. These two 
solutions do not coincide because they correspond to the 



different initial conditions. On the other hand, we may say 
that for a certain class of initial conditions, that difference 
should become irrelevant at late times i.e. the two solu- 
tions, or at least certain quantities calculated out of them, 
will coincide. But this is the same as stating the obvious 
fact that static or stationary solutions of ([TJ exist, without 
explaining whether such stationary states are indeed the 
ending point of solutions for an reasonably large class of 
initial conditions. In this light our arguments as given so 
far seem rather trivial. 

Our arguments are essentially about symmetries. The 
most convenient context to discuss them, and possibly the 
only context, is that of isotropic turbulence. We shall argue 
that the evolution of the system to the stationary phase, 
can be thought of as the gradual breaking of a larger ap- 
proximate symmetry to the smaller exact symmetry ([SJ, 
which is solely consistent with the stationary state. That 
will be realized in certain convenient cases where one may 
convince oneself that one 'watches' the system evolving as 
claimed. By the very form of ([SJ one may guess that stan- 
dard knowledge from the freely decaying turbulent flows 
could prove useful to us. 

We start by reviewing certain useful facts about the 
freely decaying isotropic turbulence. 

V. HOMOGENEOUS AND ISOTROPIC 
TURBULENCE 

A. Important quantities and formulas 

Consider homogeneous and isotropic turbulence. The 
one-direction r.m.s. value of the velocity, qi, does not 
depend on the direction, i.e. q 2 — 3q 2 . The two-point 
correlation function of the velocity is reduced to a scalar 
f(r) which depends only the distance r between the two 
points: ui(0)ui(r) = q\ f(r). The entire information of the 
two-point correlation is contained in components ui longi- 
tudinal in the direction of separation. Also the two-point 
triple correlation of the velocity can only have longitudinal 
components and it is expressed in terms of a scalar h(r) by 
ui(0)ui(0)ui(r) = qf h(r). A priori all quantities depend on 
time, for that reason time-dependence is left understood. 

Equation ([TJ with A = is the unforced Navier-Stokes 
equation describing turbulence in the freely decaying state. 
The 'Karman-Howarth equation' [l8j][lj| derived from it 
under the conditions of homogeneity and isotropy reads 

U^ = hU^ h+2 ^%)}- (18) 

In freely decaying turbulence the rate at which energy K 
is decreasing equals the dissipation rate e, expressing the 
balance of total energy in that problem. 

k = -s. (19) 

Presumably, this also holds if the viscosity v depends ex- 
plicitly on time. This fact will be useful below. 
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The integral scale, L = f Q fdr, is of the order of magni- 
tude of the dissipation length L £ . The Taylor micro-scale 
X g is defined by a differential relation involving /: 



X 2 ^l 
a Q r 2 



= -1, 



r=0 



lQvK 



(20) 



For completeness, and as we shall briefly need it later on, 
we write down the energy balance equation for the spec- 
tral densities of K and s. It is a Fourier transform of the 
Karman-Howarth equation (|18l) . see e.g. |19j |. 

(21) 



d t E(k) = -d k T(k) - 2vk 2 E{k) . 



The 'spectrum' E(k) suitably integrates to give the ki- 
netic energy and dissipation rate, K = J °° E(k)dk and 
e = 2v J °° k 2 E(k)dk. T(k) is the spectral energy flux and 
vanishes for vanishing and infinite wave numbers. Clearly 
(UHJ) follows by integrating (|2"TT) over all k, though the 
derivation of energy balance equations will be discussed 
in more detail section IVT1 



B. Scaling symmetries and power laws 

The scaling arguments given in this section are borrowed 
from [20J. The method is an application of the reasoning 
presented in section HVl 

Consider the Karman-Howarth equation (|18[) . Now per- 
form the two-parameter scaling transformation 

t -> e a t , 



e"- a q, 



r — > e r , 

q 

/-►/, 

h -> h , 



(22) 



for arbitrary constants a and b. Changing a for fixed t 
amounts to time evolution from the initial moment t. Sim- 
ilarly changing b for fixed r amounts to looking at larger 
distances. Under (|22|) equation (fl8|) becomes 



d_ 

Of 



(«?/) 
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(qlh + e^ b 2, q l d J-)}. (23) 



Consider high Reynolds numbers. Then the viscosity term 
can be dropped. We see then that the transformation (f2"2"j) 
is an approximate symmetry of (|18j) for high Reynolds num- 
bers; it can be regarded as a symmetry of the system for 
infinite Reynolds numbers. 

Consider then quantities of interest such as the kinetic 
energy K, or the integral scale L (equivalently, the dis- 
sipation length L e ). They transforms same as q 2 and r 
respectively. 

The one-parameter subgroup of the transformation ([22| 
such that 



7= ~ ) 
a 



is an arbitrary but fixed number, given explicitly by 

t-> e a t , 
r -> e la r , 

/->/, 
h — > h , 

for arbitrary a, leaves the quantities 
t^L, t 2 - 2 ^K 



(25) 



(26) 



invariant. 

Note that this way, we think of the two-parameter group 
(|22|) as one-parameter (7) family of one-parameter sub- 
groups (|23|) . Presumably, equation (j2U)) becomes identical 
to (|18l) iff a — 2b = that is a — 2-fa — 0. This means that 
the subgroup 7 = | is an exact symmetry of the freely 
decaying turbulence. In other words, the larger symme- 
try (|2"2"|) for infinite Reynolds number breaks down to its 
subgroup 7 = i for finite Reynolds numbers. 

Each symmetry (|25[) is essentially time-evolution. Fol- 
lowing the arguments of section IIVI we conclude that at 
adequately late times 



L = constant t 1 



K = constant t 21 



(27) 



Thus we have obtained certain power laws for the late be- 
havior of the length scale and kinetic energy in freely decay- 
ing turbulence. The law for the dissipation rate e follows 
immediately by (|T9|). 



£ = constant 



(28) 



Then, by (ITfJ)) the law for dissipation length L e turns out 
consistent with that of L, as course it should. The law for 
Re^ also follows from (fT6j) : 



Re^ = constant t 



27-1 



(29) 



Summarizing, each value of 7 defines a subgroup of the 
full symmetry group (f2"2"]l for high Reynolds. Given a 7 
the time-dependence of various quantities takes the form 
of specific power laws. A priori not fixed without addi- 
tional conditions, the exponent 7 may be given an addi- 
tional physical interpretation. Assume that for low wave- 
numbers k the spectrum E{k) is of the form 



E(k) = Ck a + o(k a ) . 



(30) 



for some constants C and a. Given the dimension of E(k) 
and k and the constancy of C, this relation is invariant 
under {25) iff 

7 = -^- (31) 



That is, the subgroup ([25)) is fixed by the small wave- 

(24) 

number behavior of the spectrum of the specific class of 



7 



flows. It may be argued, see e.g. Ref. [2l|, that C is ac- 
tually constant as long as 1 < a < 4; also the case a = 4 
holds marginally. That is, in those cases C is fixed by the 
initial conditions. 

Decay exponents are usually expressed in terms of n = 
2 — 2"/, which is the kinetic energy decay exponent, K ~ 
t~ n . About the value of n there are well known sugges- 
tions. They depend upon the identification of C with 
quantities which are conserved under certain conditions. 
Kolmogorov [22j and Batchelor f23|. based on the conser- 
vation of the Loitsyanky integral [24|, derived 7 = j i.e. 
n = Saffman [25| set forth the hypothesis that the vor- 
ticity and not the velocity correlator is an analytic func- 
tion in spectral space, by which rediscovered the a = 2 
spectrum and Birkhoff's integral (26j and derived 7 = | 
i.e. n = |. Experimentally [27]][28| both values of the de- 
cay exponent n are acceptable. The value n = 1 has also 
been suggested by other theoretical considerations, for high 
but finite Reynolds numbers [2!| and as the limiting value 
of the decay exponent for infinitely high Reynolds num- 
bers [3(il - [32l ]; this solution first appeared in (33j . 

The n = 1 decay solution for finite Reynolds numbers 
of Ref. [2^] can be obtained by recalling an observation 
given above, that for finite Reynolds numbers the symme- 
try (|2"2"j) , essentially associated with infinite Reynolds num- 
bers, breaks down to its 7 = h subgroup at finite Reynolds 
numbers. That means n — 1. Presumably, by (|29|) . the 
Reynolds number is constant for this solution. 

The n = 1 decay law may also be obtained in another 
way, which gives the chance to make an additional com- 
ment on the analysis presented in this section. The Tay- 
lor microscale X g transforms as a length, same as r, ac- 
cording to the equation on the left in (|20|) . That means 
that X g — constant t 1 . On the other hand the equation 
on the right in (J2DJ and the laws <[27j) and (|28j imply 
\ g = constants. The reason why there is no discrepancy 
is because regarding L as finite and Reynolds number vir- 
tually infinite, for the symmetry (f2"2"j) to hold, means that 
X g is virtually zero. Put differently, if we want to think 
of the previous analysis as applying also to high but finite 
Reynolds numbers, then we must restrict ourselves to scales 
much larger than Taylor microscale. It is then no accident 
that the power laws (|27|) can also be produced by models 
deriving from self-similarity of turbulence with respect to 
the integral scale L, as we shall discuss in section fVIII On 
the other hand if we want finite Reynolds and to take into 
account scales of 0(X g ) or less, then it must be 7 = \ i.e. 
n = 1. 

Being such a direct implication of the arguments in this 
section, one may wonder why the n = 1 decay solution is 
not observed experimentally even for the highest Reynolds 
numbers (equivalently, as 7 = | means a = 1, a small 
wave-number spectrum E(k) oc k has not been verified). 
The arguments possibly fail on the very part where one 
expects independence from the initial conditions. That ex- 
pectation might be in a better shape the higher, still finite, 
the Reynolds number is. This is why in the best case the 
n = 1 solution can possibly be regarded as describing well 



decaying turbulence for very high Reynolds numbers. 

In the next two subsections we come to the problem of 
interest. The discussion parallels in some sense our previ- 
ous remarks: Going from the infinitely high to any lower 
Reynolds number the larger symmetry (|22[) breaks in this 
case down to its exact subgroup 7 = associated with 
the linearly forced turbulence, the exact symmetry Q we 
started our discussion with. But, unlike the freely decaying 
case, in our problem a large length scale and a Reynolds 
number scale are necessarily present, eventually forcing the 
system towards the 7 = evolution. That amounts to 
reaching the stationary state. 



C. Linearly forced isotropic turbulence 

Consider linearly forced turbulence in the description 
given by equation (|5|), which let us state again: 



vV. 



The analogue of the transformed Karman-Howarth equa- 
tion (1231) for late times t' 3> A^ 1 reads now 



dt 



Mf) 



1 d 



l fn' 3 h' + p- 2 ^ a — 9n' 2 9 f "0 

[q.h+e M 2q x ^Jj 



We find of course again that the subgroup 7 = of the 
group (|22t i.e. the group ([9]), is an exact late-time sym- 
metry of the linearly forced turbulence. For very high 
Reynolds numbers the larger group (|22|) is a good approx- 
imate symmetry of the system. The evolution laws for the 
dimensionful quantities are necessarily similar to those of 
the freely decaying case, 



V = constant t' 1 , 
K' = constant t' 27 " 2 . 
e' = constant t' 2l ~ 3 . 



(32) 



while the dimensionless Reynolds number has a different 
power law due to the time- varying u 1 : 



Re'r = constant t' 21 . 



(33) 



Consider a flow that starts off with velocities of order uq 
and a box of size I such that uq 3> Al. Equivalently the 
turn-over time is much smaller than forcing time scale A^ 1 , 
that is, l/uo <C A^ 1 . 

Given A, I and v there is a naturally defined Reynolds 
number in the problem: 



Al 2 

Re a = 

v 



(34) 



That is, the condition l/u -C A^ 1 can be rephrased as 
that the flow starts off with a very high Reynolds number, 
Rei > Re^- 

Consider then times t' such that 1/uq <C t' <C A^ 1 . Look- 
ing at the previous equation we understand that for those 
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times the turbulent flow is merely freely decaying with con- 
stant viscosity v. If all previous inequalities hold strongly 
enough, then there will be time for the flow to evolve ade- 
quately towards its developed stage. That means that the 
quantities describing turbulence will evolve according to 
the power laws (|32t and (|33|) . 

When t' becomes of order of A -1 'linear forcing' kicks 
in. Now one should recall that Reynolds number is always 
decreasing, therefore some time before or after that mo- 
ment it will drop enough so that the viscosity term cannot 
be neglected. That means that from that moment on the 
group (|22[) is not much of a symmetry anymore: the only 
symmetry remaining is its subgroup (0) corresponding to 
7 = 0, which is exact and therefore holds at all times. Vis- 
cosity now decreases with time therefore energy will be dis- 
sipated with an ever decreasing rate. We may then picture, 
very roughly, the flow evolving by going through stages of 
smaller exponents 7, following simple of less simple laws 
parameterized by it, eventually reaching the specific value 
for which a subgroup (j25p is an exact symmetry of the sys- 
tem: 7 = 0. 

Let us recapitulate. There is a symmetry existing in 
the system for high Reynolds numbers. A priori allows 
for arbitrary values of the parameter 7. This symmetry 
breaks down to its subgroup 7 = when Reynolds drops 
enough. This is an exact symmetry of the system, therefore 
holds always. The final state must be the one respecting 
that exact symmetry. The power laws (|32p and (1331) imply 
that L' and Re^ are constant, e' = constant <'~ 3 and K' = 
constant i' -2 . The transformations (jT^J) and (|T71) back to 
the original variables of equation ([1} show that everything, 
that is K, e, L e and Rez,, is constant. We have reached 
the stationary state. K' and e', which are time-dependent 
in the description © of the problem, are related by the 
equation (fT9|) written down for the primed quantities: 



The l|32p and (f3"3")) power laws for 7 = and the transfor- 
mations (H3J) translate that relation to 

2AK = e . (36) 

That is, at the stationary state energy production balances 
exactly dissipation. 

If initially Reynolds number is not very high, there will 
be a strong mixing of the phases of freely decaying and 
'linear forced' turbulence. The transition between them 
is then too complicated to describe. Though one cannot 
argue against the possibility the system having an entirely 
different behavior, it seems reasonable that the same mech- 
anisms which lead the system to the power laws consistent 
with 7 = will do the same thing, in far more complicated 
way. 



D. Effects of the finite domain 

A very amusing thing we observe is that, by pip , the 
value 7 = corresponds to a = 00. This indicates that 
the power law behavior of the spectrum E(k) for small 
wave-numbers degenerates, and should be replaced by some 
other, much faster decreasing law, perhaps some kind of 
exponential. 

There is a good reason why we expect that. The sys- 
tem ((T|), or ©, is solved in a domain of some finite size 
I. An infinite size I is meaningless: In the absence of 
another large length scale, this means that the total en- 
ergy production rate in the domain depends on it and di- 
verges 1 , 2AKpl 3 — > 00. Also large I essentially means a 
large Al compared with any specific initial condition uq: 
an infinitely large I is equivalent to initial conditions uq 
infinitely close to zero. Now finite size means that there 
are no wave-numbers k between 0{l~ l ) and zero, therefore 
any continuous approximation of the spectrum E(k) must 
fall very rapidly for kl smaller than 0(1). 

There is a major implication following the presence of a 
finite size domain. Its fixed size I breaks the symmetry i22\ ), 
as the presence of a fixed length says that it must be b = 0. 
That is, the domain size breaks the larger symmetry (|2"2"j) 
down to its subgroup 7 = 0, the exact symmetry. 

We may now think of the evolution of the flow from an- 
other point of view, that of the integral scale. As long as 
the integral scale L' is small compared to I the group (f2"2"|) 
is a fairly good approximate symmetry. Then L' increases 
with time as ~ t' 7 . As L' grows larger, (|2"2"|) is a less and 
less good approximate symmetry. As before, we may then 
roughly picture the flow as going through stages of smaller 
7 reaching the stage with 7 = which is consistent with the 
exact symmetry. This means that V will become constant. 
The natural order for that constant V (as well as L, recall- 
ing that L = L' by the transformation (|17p) is the domain 
size I. As mentioned in the Introduction, DNS have shown 
that specifically L E = I within a few percent error 

We may note that for high Reynolds numbers the viscos- 
ity term can be anyway neglected. Therefore the (J3U) and 
(I33p power laws for 7 = and the arguments given above, 
make sense for very high Reynolds numbers for freely de- 
caying turbulent flows in a finite domain. The difference in 
the 'linearly forced' case <JSj> is that those power laws are 
associated with an exact symmetry of the equations and 
of course hold also for lower Reynolds numbers. In con- 
trast, as we have seen in section IV Bl the exact symmetry 
in the freely decaying turbulence requires 1 — \- Linearly 
forced turbulence behaves in a simpler way than the freely 
decaying, in this sense. 

One should note that by the time linear forcing kicks 



The density p is regarded as constant, which for infinite I means 
also infinite mass p/ 3 . But that it is not a priori harmful, because 
what matters is the largest scales where turbulent motions correlate 
and the total energy produced in the system. 
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in, the time-dependence of viscosity v' changes the form of 
the Reynolds number power law. During times that linear 
forcing is already at work but the Reynolds number is still 
very high compared to Re a and decreasing, the Reynolds 
number is given by the power law (|3"5|) and decreases due to 
a decreasing 7 as we argued in the previous section. That 
is, when 7 eventually vanishes the peculiar 'decaying' tur- 
bulent flow ([5]) reaches a peculiar kind of stationarity: its 
Reynolds number becomes constant i.e. turbulence as such 
is not decaying at all. This reflects of course the association 
of this state to an scaling exact symmetry of the full viscous 
equations. Recalling transformation (1171) . Re L = Rcl, this 
is also the Reynolds number of the linearly forced turbu- 
lence described in the variables of (jTJ) . 

That final value of the Reynolds number can be easily 
estimated by the existence of a natural scale Re^, defined 
in (f34|) . At the stationary state energy production balances 
dissipation, 2AK — e, which implies that Re^ = AL 2 E /(Av) 
i.e. Re^ indeed sets the scale for the Reynolds number at 
the stationary state. In fact, as mentioned in the Intro- 
duction, DNS show that Re^ = -jRe^ within a few percent 
error \Tl\ . 

We conclude that the finiteness of the domain emerges 
as a crucial factor in understanding the flow evolving to a 
stationary state. Now the domain, we often referred to as 
the 'box', should be understood more carefully in terms of 
periodicity. This is what we discuss next. 



VI. THE STATE OF ISOTROPY 

In the previous section we restricted ourselves to flows 
obeying the conditions of homogeneity and isotropy. Let 
us review what that involves. Contract (TTJ) with u and 
average. After a little re-arranging we have 



dK 
~dt 



= -c + 2AK 



-vSJ K - V ■ J . 



(37) 



P 



Homogeneity alone makes all locally defined correlators, 
such as K or Ji, independent of the position in space. That 
is, the last two terms in the last equation vanish. Of course 
everywhere isotropy means homogeneity, so if the flow is 
assumed isotropic those two terms go again away. We are 
then left with 



K = -e + 2AK . 



(38) 



In a stationary state we get e — 2AK, whose form we 
already anticipated on dimensional grounds deriving the 
final value (|34p of the Reynolds number in terms of the 
box size I. 

We should now recall that the box is a cubic domain to- 
gether with periodic boundary conditions we impose on all 
fields. The periodic boundary conditions can be given some 
enlightening interpretations. One way to think of them is 



that we solve the Navier-Stokes in an infinite medium im- 
posing periodicity I on the field u(x, y, z) in all three direc- 
tions x, y and z. That in turn means that homogeneity is 
not a priori broken: a box such that u(x, y, z) satisfies pe- 
riodic boundary conditions can be drawn anywhere in the 
infinite medium. By the periodicity we apparently restrict 
ourselves to special kinds of flow such that there is an up- 
per bound to the size of eddies or any spatially periodic 
structure in it. 

Another way to think of the periodic boundary condi- 
tions is to interpret them as mere single- valued-ness of the 
fields while identifying the points of the boundary of the 
cube where the fields are supposed to be equal. This may be 
pictured if we go one dimension down. If we take a square 
and identify the opposite sides we get a topological torus. 
A torus is a perfectly homogeneous space without bound- 
ary which is non-trivial globally and it is not isotropic. 
Imposing periodic boundary conditions on the cube means 
we essentially solve Navier-Stokes equations on the three- 
dimensional analogue of such a space, the three-torus. 

An implication of periodicity and its peculiar nature is 
that an analogue of the equation holds without as- 
suming point-wise homogeneity of the flow. Let us denote 
by (X) the spatial average over the volume of the box of 
a quantity 2 with an ensemble average X. Averaging that 
way we get an equation similar to Q37[) for box-averages 



d(K) 
dt 



= - £ 



Vb, 



2A(K) + 

— [ {uVK-J}-da. (39) 

30X Jbdy 



The last term is a surface integral over the boundary 
surface of the box. This integral is a sum of 



(vd x K - J x )dydz 



x—l 



[vd x K - J x )dydz (40) 



x=0 



plus two more analogous pairs of terms for the other two 
directions. As all fields are manufactured out of correla- 
tions of the field u which satisfies periodicity u(x, y, z) = 
u(x + l,y,z), and the same for the other two directions, 
any pair of terms such as (|40| vanishes exactly and identi- 
cally. In other words even if point- wise homogeneity is not 
assumed, periodicity implies that 



d(K) 
dt 



■(e) + 2A(K) 



(41) 



must hold exactly. 

Alternatively, the vanishing of the boundary term in 
(|3"9")l follows automatically under the interpretation that 
we solve our problem on a three-torus: there simply is no 
boundary, x = and x = I describe the same surface some- 
where on the three-torus. Whatever the interpretation of 



That is, for example, (K) means simply (~U'll), not (|u-u). This 
is a little confusing but allows for a more compact notation. 



10 



0.50 




0.20 -I 1 1 1 1 i 1 

20 40 60 80 100 120 140 

3 At 



FIG. 2: Time evolution of the diagonal components of the nor- 
malized Reynolds stresses, (u a u a )/q 2 (no sum), from the same 
run producing the time-series in fig. [T] 

the boundary conditions, one ends up with (|4ip without 
assuming homogeneity of the flow. 

This is a good thing to know. The box-averaged correla- 
tors (X) are the actual observables in the DNS. Now, what 
is their relation to the ensemble averages XI Assuming 
that X are meaningful and under an ergodic hypothesis, 
X are represented by averaging (X) over suitable and ade- 
quately large intervals of time. That first of all means that 
although the motion of (X) is not the same as that of X 
it should nonetheless be bounded and appear as fluctuat- 
ing around hypothetical stationary values. This is what it 
is observed, fig. Q] Those values should of course be the 
values of the correlators X. 

Contemplating the far more complicated motion of the 
correlators (X) compared to that of X, one quickly realizes 
that both kinds of correlators obey the same basic dynam- 
ical equations. Thus their difference lies somewhere else. 
The simplicity of the motion of X follows from their in- 
dependence from the initial conditions of the flow and the 
symmetries of the dynamical equations and of the domain. 
Correlators (X) are not a priori independent of the initial 
conditions of the flow and none of the arguments of section 
IfVI applies to them. Therefore there are many more and 
more complicated solutions (X) than X. 

In particular, in section IIVI it was emphasized that the 
symmetries of the cube, combined with the symmetry 
force a fair amount of isotropy on the solutions X. That 
will hold only on the average for the correlators (X) . This 
is the actual result of numerical simulations, fig. [2J 

We consider the fluctuations of measures of isotropy 
important for the following reason. The 'fair amount of 
isotropy' exhibited by X is not so harmless itself: Rea- 
sonably, any anisotropy is inherited by the correlators (X) 
and it is enhanced in their description of turbulence. The 
origin of any enhancement of anisotropy is that in linear 
forcing there is no intrinsic large scale at which we feed the 
flow energy in some isotropic manner; the large scales are 
set by the domain itself and the scales comparable to its 
size are necessarily not isotropic. Subsequently, anisotropy 
is produced and maintained at those and smaller scales 



through forcing and cascade. Although the motion of the 
correlators X suggests that the system cannot be kicked off 
balance, the produced anisotropy will cause relatively large 
fluctuations of all quantities (X) describing turbulence. 

In the next section we will try to lend some quantita- 
tive support to this intuitive picture. If deviations from 
isotropy are related to the fluctuations of all quantities then 
any fluctuations should vanish in a perfectly isotropic set- 
ting. That may give us a sense of what happens when 
those fluctuations are continually generated. As we shall 
see the analysis is interesting in its own right as it reveals 
rather important dynamical properties of linearly forced 
turbulence. 



VII. SELF-PRESERVING TURBULENCE AND 
STABILITY OF STATIONARITY 

Studying the fluctuations around the stationary state is 
equivalent to studying the stability of that state as a fixed 
point of solutions, in the statistical sense. In fact, this is 
an alternative way to look at the main problem we have 
been concerned with, the stationary state as an attractor 
of solutions. Of course, such an analysis is a very difficult 
thing to do unless we resort to some suitable simplification. 

According to the plan set at the end of the previous sec- 
tion, we shall assume that the flow is isotropic. Therefore 
all correlators involved, which can be thought of either as 
ensemble correlators X or box-averages (X), are assumed 
to have the properties required by that condition. We may 
then investigate the fate of any deviations away from the 
stationary state if the flow evolves remaining isotropic. 

The evolution laws derived in the previous sections can 
be alternatively derived in isotropic turbulence from models 
relating the kinetic energy K and dissipation rate e. These 
models can be deduced on dimensional grounds, or more 
systematically by self-similarity arguments, which are fairly 
equivalent to the scaling arguments given here. The latter 
date back to the work of von Karman and Howarth [34j and 
Batchelor j23[, see also [35| 36]. One looks for self-similar 
solutions of the equations w.r.t. a single length scale L(t), 
'self-preserving' turbulent flows. Assuming that the larger 
scales of the flow are evolving in such a self-preserving man- 
ner, one chooses L(t) to be the integral scale and one ob- 
tains a closed system of equations for the variables K and 
s. That simple model can also be regarded as describing 
self-preserving turbulence of all scales but for infinitely high 
Reynolds numbers, essentially for inviscid flow. 

One can straightforwardly apply the same arguments in 
the linearly forced turbulence. The spectral energy balance 
equation (l2TT) becomes 

d t E(k) = -d k T(k) - 2vk 2 E(k) + 2AE(k) . (42) 

The origin of the additional term should be clear. Then 
the self-preserving development of the larger scales of the 
flow implies, via standard steps which can be found e.g. in 
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35] [36J , the model equation 



de 
dt 



-Ci 



K 



■ c\Ae , 



(43) 



where c\ — 3 and Cf is a dimensionless constant. Apart 
from the value of ci, this equation could also have been 
guessed on dimensional grounds upon requiring its r.h.s. 
to be built out of e and K and A and be linear in A. 

Integrating (|4"2")l over all wave numbers we obtain again 
the exact equation (|38[) . which we write down again for 
convenience, 



d JL = -e + 2AK. 
dt 



(44) 



The system of equations (14"4"|) and (|4"3")l is consistent 
with a static solution only for 2Cf — c\. The special 
case = c\/2 = 3/2, predicted by large scale self- 
preservation, implies that L £ = constant at all times the 
model holds. This is consistent with the general idea about 
it. The large scale self-preservation model equation (143)) 
and the value =3/2 will emerge again from a different 
perspective in section IVIIII The model can be easily solved 
exactly and indeed predicts that the flow approaches sta- 
tionarity exponentially fast for all > 1 (the case = 1 
is trivially consistent with stationarity) . 

A more elaborate analysis of the evolution of isotropic 
turbulence has been presented in the past in the 
Refs. [13 dllHII] [11. In those works the self-similarity hy- 
pothesis is applied at the viscous equations of the flow i.e. 
self-preservation is required to be true for all scales of tur- 
bulence for finite Reynolds. In the terminology of Ref. [29], 
self-preservation is complete. An implication of this re- 
quirement is that the self-similarity scale is the Taylor mi- 
croscale X g . 

From the point of view of the linearly forced turbulence 
all that sound very relevant and interesting for the following 
reasons. First, the linearly forced turbulence comes to the 
intelligible part of its course when its Reynolds number 
approaches the value (13"4"|) which need not be very high 
at all; second, energy is generated uniformly at all points 
in the domain and it feels that all scales play a role in 
approaching or maintaining stationarity; and third, in this 
problem there is a natural scale for the Taylor length \ g . It 
is the scale at which energy production balances dissipation 
in spectral space, as can be seen by equation (|4"2"j) : 




(45) 



This is designated as a Taylor microscale because the sta- 
tionary state value of the Taylor microscale, A gs , is of that 
order: 



A„ 



V5A, 



(46) 



This follows from the definition (|2U| of X g and the station- 
ary state total balance of energy production balances dissi- 
pation, 2AK = e. For these reasons the Taylor microscale 



may be regarded as playing a particularly significant role in 
the dynamical aspects of linear forcing, perhaps quite more 
significant than in the freely decaying case. [On the other 
hand, as everything turns out approaching constancy, even- 
tually the integral scale might be used as a self-similarity 
scale, a choice associated with the model (|43|) . providing 
a more crude and late-time description of the evolution of 
the system.] In any case this choice does provide a closed 
two-equation model with some interesting properties. 

We may then proceed as follows. There is another exact 
equation which we may use along with (|44[) . One way to 
derive it is to start from the Karman-Howarth equation for 
linearly forced isotropic turbulence 



d_ 

dt 



(«?/) 



r 4 dr I 



(q\h 



2vq\ 



0A 
dr J 



}+2Aqff, (47) 



applying the definitions (l4"9")l below. (The spectral energy 
balance equation (|42p is a Fourier transform of (|T5|) .) Al- 
ternatively, and more instructively, we can do everything 
from scratch by differentiating s w.r.t. to time using its 
very definition as an ensemble or box-average correlator. 
Then, employing the Navier-Stokes equation (|T|) and ap- 
plying the condition of isotropy on any arising correlator 
one arrives at 



ds 
dt 



7\S\ 



3/2 _ 7<?£ 

TED 15K 



■2As, 



(48) 



where S (the velocity gradient distribution skewness) and 
G are defined by 



G - Xg 



9V 
dr 4 



(49) 



where / and h are the two-point double and triple point 
correlations of the velocity defined in section lV Al Equation 
(|48|) can also be derived by multiplying (l42t by 2vk 2 and 
use formulas equivalent to (|49j) and (|20j) in wave number 
space. 

The system of equations (|44l) and (|4"5)l is not closed, the 
dependence of S and G on K and e is unknown. Assume 
now that at some moment to the flow becomes self-similar 
with a (time-dependent) similarity scale A . That means 
/ and g are functions of the dimensionless coordinate r/Ao 
alone, modulo a possible dependence on the initial condi- 
tions at to. Now ([21?]) tells us that A /A ff must be a con- 
stant, depending only on the initial conditions at to- Thus 
the similarity scale is indeed the Taylor microscale. Then 
by (|49|) we have that S and G are constant and equal to 
the values they have at to'- S = So and G = Go. Now the 
system (|4"4"| and (143)) is closed and we may study it. 

Let us denote the stationary state values of the dissipa- 
tion rate and kinetic energy by e s and K s . Of course they 
are related by e s = 2AK S . We would like to study the 
stability properties of e — e s and K — K s as a complete 
self-preserving solution of the system of equations (|44]) and 



It will be convenient to define the quantity 

JGo 
15 
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(50) 
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First of all, equation 



implies that 



61 3At 81 



a/2 



2A- 



which implies that 



7|5 | 



9 >1 



(.9 - 1) , 



(51) 



It is useful to relate the value of g to the Taylor-scale 
Reynolds number Re a = (-j^Rei,)*. By (JT5J) we find that 
its stationary value Re^s reads 



ReAs = 



30 
7iSo 



(9-1) 



(52) 



Define now small fluctuations £ and £ of e and K around 
their stationary values: 



e = e s (1 + , A" = JT, (1 + C) . 



(53) 



Inserting these expressions into (|44|) and (|48|) and keeping 
only linear terms we obtain the following system of equa- 
tions: 



^ t =-A{l+g)i + 2AgC,. 
^ = -2A£ + 2A(. 



(54) 



Its eigenvalues T read 



1 



-A 



(-(3-l)± V(.9-l)(.9-9)) 



(55) 



By g > 1 we see that the real part of both eigenvalues is 
always negative. Fluctuations around the stationary state 
die out exponentially fast. That is, modulo finite domain 
effects, the stationary state is stable as a complete self- 
preserving isotropic solution. We may also view this result 
as providing further evidence that the stationary state is 
the natural final state of the linearly forced turbulence. 
[Presumably, one may observe that exponentially fast ap- 
proach to the stationary state is also the prediction of the 
simpler model (|43[).] 

The previous analysis can be alternatively understood as 
follows. In order to derive the previous results we have as- 
sumed perfect isotropy. A reasonable assumption about the 
deviations from isotropy is that they originate from scales 
of order That means, according to our conclusions in 
the previous section, that the same can be said about the 
fluctuations around the stationary state. That is, one may 
attribute the generation of fluctuations to the interaction 
of the larger eddies with the periodicity i.e. the restriction 
to their size. Then, through both forcing and cascade, fluc- 
tuations are generated at all scales from I down to a certain 
scale where isotropy becomes a good approximation. There 
things are different. We may define correlators as spatial 
averages (X)y over volumes V smaller than that maxi- 
mum isotropic scale i.e. within these volumes turbulence 
is isotropic (meaning homogeneity as well) to a good ap- 
proximation. Then K and e understood as spatial averages 




3At 



FIG. 3: The time evolution of the dissipation rate £ which is 
shown in fig. [1] is shifted in this figure by one unit of dimen- 
sionless time 3At. 

(X)v obey similar equations to those studied above. The 
entire previous analysis goes through. That finally means 
that at adequately small scales the fluctuations are strongly 
suppressed, but at all higher scales are maintained through 
forcing and cascade. The maximum isotropic scale should 
be (very) roughly related to the characteristic Taylor mi- 
croscale of linear forcing A ,4 = {v/A)?, as below that scale 
dissipation becomes stronger to energy production. 

We may investigate the linear system ((54)) a bit further. 
Though this system meant to serve us mainly for qualita- 
tive considerations, regarding the stability of the constant 
solution K = K s and e = e s , there are some amusing re- 
marks to be made about it solutions on the quantitative 
side. In the range 1 < g < 9 the eigenvalues T are complex 
numbers. If we take for definiteness \Sq\ = 0.5, that means 
that when Re\ < 69 the fluctuations are damped oscilla- 
tions. [Presumably, this emergence of oscillations is a qual- 
itative difference between the complete self-preservation 
model and the simpler model (|43l) .] Inserting the solutions 
C = Coe r 'and £ = £oe rt , for positive frequency into any of 
the equations (fM)) we obtain the phase difference and the 
relative amplitude of e and K: 



where <fi is given by 



tan < 



y/(g-l)(9-g) 



(56) 



(57) 



As expected the dissipation e evolves with a phase delay 
w.r.t. the kinetic energy K and the energy production 
2 AK . This corresponds to a time-delay (/>/|Imr|. The pe- 
riod of these damped oscillations is of course 27r/|ImT|. 

In fig. [1] we plotted the energy production 2AK and dis- 
sipation rate e against time in units of (3^4) _1 . Now let us 
shift the evolution of dissipation by one unit of time to off- 
set its delay. The result is given in the fig. [31 One observes 
that after that shift the complicated oscillations appear in 
phase to a considerable degree of accuracy. Curiously, the 
time-delay 0/|Imr| in units of (3A) _1 decreases from 1.5 
to 0.5 in the range 1 < g < 9. Also the period 27r/|Imr| is 
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15 
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(c) 



FIG. 4: Time evolution of the kinetic energy (solid line) and 
dissipation (dashed line) normalized by their stationary state 
values from numerical solutions of the system (|58fl . Figure (a) : 
initial conditions K(0) = 1.3 and e(0) = 1.2 for g = 2.2 that is 
Re As ~ 10. Figure (b): K(0) = 0.01 and e(Q) = 0.01 for the 
same value of g. Figure (c): K(0) = 0.01 and e(0) = 0.01 for 
g = 9.9 that is ReA S ~ 85. 

roughly an order of magnitude higher than (3A) _1 for most 
values of g, which as a number is not in disagreement with 
the picture in fig. [3] Given that these numbers derive from 
a model which does not interact with the source of the fluc- 
tuations, it seems interesting that the oscillations it implies 
may encapsulate certain features of the actual fluctuation. 
There certainly is no identification between the actual fluc- 
tuations and those oscillations. For example, when g ~ 9 
that is Re\ s ~ 70 the damped oscillations are replaced by 
a purely decaying exponential, a qualitative change in the 
behavior which cannot be traced in the DNS results of the 

Refs. [Ill [ni [Hi- 

The previous remarks derive from the quantitative char- 



acteristics of small fluctuations, and we may have over- 
extended the applicability of the related formulas. Arbi- 
trary fluctuations arc described by the solutions of the full 
non-linear model ([35]) and f4"8"j) . This needs to be solved 
numerically. In terms of the dimensionless (hatted) kinetic 
energy, dissipation rate and time, defined respectively by 
K = K s K, e = e s i and t = 2At, the non-linear model 
reads 



dK 



(58) 



The parameter g is related to ReAs by (|52p and we again 
take for definiteness |So| = 0.5. 

The system ([55)) is solved using the software Mathemat- 
ica. We consider a few specific cases. First, the difference 
of the initial conditions from the stationary state values is 
such that to imitate the size of the observed fluctuations. 
This is shown in fig. Second, the kinetic energy K and 
dissipation rate e start off from very close to zero, shown in 
fig. Sb. For those two cases we have chosen an adequately 
small Reynolds number so that oscillations to be visible. 
Finally we consider the effect of higher Reynolds numbers. 
An evolution of K and e for Ke\ s ~ 85 is shown in fig. [4b. 

The result is that the picture is not qualitatively different 
than the one obtained from the small fluctuations. In the 
figs. 0k andllb, the time-delay of the dissipation relatively 
to the kinetic energy and the period of the damped appear 
essentially as predicted previously, and the oscillations of 
the dissipation are consistently larger as implied by equa- 
tion (|S"f)|) . On the other hand fig. [3b> shows a particular 
behavior of the non-linear solutions: if the initial condition 
is far away from the stationary state values the system un- 
dergoes large fluctuations before settling to those values. 
The fig. 0b shows that increasing the Reynolds number any 
wiggling of the curves due to oscillatory behavior dimin- 
ishes to extinction, which is again what we expected. 



VIII. DISCUSSION 

Direct numerical simulations of turbulence forced by the 
linear forcing scheme exhibit a not entirely expected sta- 
tionary late-time state. The stationary phase is essentially 
quasi-stationary: all quantities have relatively large fluctu- 
ations, though their time-average can be predicted fairly 
well. In the present work we have attempted to under- 
stand how these phenomena are rooted in the properties of 
the system. This was done by using the symmetries of the 
dynamical equations of the problem, as well as the sym- 
metries associated with the boundary conditions i.e. the 
size and the symmetries of the cubic domain; also, using 
special dynamical properties of the system derived under 
usually employed conditions such as exact isotropy or self- 
similarity. In this problem there are few and specific scales: 
the domain size I, the constant rate A and the viscosity v. 
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Out of them derive a Taylor microscale \a and a Reynolds 
number Re^- These quantities control the major (intelligi- 
ble) features of linearly force turbulence evolution. 

The importance of the finiteness of the domain and its 
effects cannot be over-emphasized in the linearly forced tur- 
bulence. In a limited bandwidth forcing scheme, determin- 
istic or stochastic, the inverse wave numbers at which one 
forces the flow imitate, very roughly, the scale of a phys- 
ical stirring of an incompressible fluid existing in slightly 
larger 'box'. In linear forcing there is no such intrinsic 
scale. This simplifies things in some sense because there is 
no interaction between the forcing and domain size scales. 
On the other hand it is left entirely to the domain to set 
the large scales, becoming an essential part of the forcing 
itself. Also the large scale is introduced geometrically as 
a matter of size and not dynamically as in the bandwidth 
schemes, and there is no actual control over the extent 
forcing is consistent with isotropy. Turbulence is expected 
to behave quite differently under linear forcing than under 
a bandwidth scheme. There some additional interesting 
properties of linear forcing we have not yet commented on. 
These properties can be associated with the effects of the 
finite domain size, and also show an at least formal affinity 
of the linear forcing to freely decaying turbulence, than to 
the bandwidth forcing schemes. 

Denote by Aui the longitudinal velocity difference. The 
second and third order structure functions are related to 
the correlation functions / and h, introduced in section 
El by (Aui) 2 = 2qf{l - /) and (Auj) 3 = 6qfh. For ade- 
quately high Reynolds nu mbers th ere is a range of distances 
(the inertial range) where (Au/) 2 = C^er) 2 / 3 , where C2 a 
constant. Consider first decaying turbulence. It evolves 
according to the power laws (|27)) . the integral scale is pro- 
portional to i 7 . The law for e can be deduced. It is then 
straightforward to show that they satisfy the K — e model 
equation ([4*3"]) for 



replacing 



3-27 
2-27 ' 



(59) 



and of course A = 0. Using the Karman-Howarth equa- 
tion (TT5)) it then straightforward to show [39] 4Cj that for 
very high but finite Reynolds numbers, and within the in- 
ertial range (more specifically as long as r/X g is a number 
of O(l)), the two-thirds law of the second order structure 
function implies specific finite Reynolds number corrections 
to the four-fifths law of the third order structure function, 



of <3(Rc 



-2/3\ 



(A Ul ) 3 = 
5 x 153 



The result is 391 4C 



-erx 



(60) 
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Consider the same question in the linearly forced tur- 
bulence. One may follow the same steps, starting from 
the Karman-Howarth equation with linear forcing, equa- 
tion (|47p. One finds a result entirely similar to (|60[) upon 



Kde 
e^dt 



3AK 



(61) 



Observe now that if we think of the r.h.s. of this sub- 
stitution as a constant, then we re-discover the model 
equation f|43|) : the constant is what we denoted there by 
C^. Equation (|4"3"]) is derived assuming self-similarity (self- 
preservation) of the larger scales of turbulence with re- 
spect to the integral scale L for high Reynolds numbers, 
in both the linearly forced (A ^ 0) and freely decaying 
case (A = 0). In all, by self-preservation we obtain a sim- 
ilar result of the form ((60)) in both kinds of turbulence, 
differing only in the value of the constants and C e . On 
the linearly forced side, self-preservation requires Cf = 3/2 
and equation (|4"3"]) and (|4"4"]) require that L = constant. At 
first sight there is no such restriction on the freely decaying 
side. In all, there appears to be a correspondence between 
linearly forced and freely decaying turbulence, though this 
correspondence appears inexact. 

Now if we require — C e then by ([591 we have that 
7 = 0. In other words, if the decaying turbulence evolves 
according to L ~ constant (and K ~ t~ 2 ) then its struc- 
ture function expression ([60]) is exactly similar to that of 
the linearly forced turbulence. That is, the correspondence 
between the two flows can be exact. 

The K ~ t~ 2 evolution is too fast compared to the usu- 
ally observed decay laws, discussed in section IV Bl Such 
power laws can be reproduced if choose the constant C e to 
be different than 3/2, a fact regarded as an imperfection 
of the correspondence in the Ref. [lj|, where it was first 
pointed out. On the other hand the origin and the nature 
of the correspondence seem to have been overlooked in [15j . 

The key role is played again by the finiteness of the do- 
main. As emphasized in section IV Dl a container is a nec- 
essary thing when turbulence is linearly forced. Lacking 
an intrinsic length scale, linear forcing essentially requires 
a large scale to be provided by the boundary conditions. 
It is therefore not much of a surprise that similarities be- 
tween linearly forced and freely decaying turbulence are 
more detailed when the decaying side evolves in a way con- 
sistent with the existence of a container: For adequately 
high Reynolds numbers that means L ~ constant (and the 
rest of the (g7|)-([2S|) power laws for 7 = 0). Then the math- 
ematics of self-similarity of turbulence with respect to the 
scale L imply exactly the same formula (|60| for both kinds 
of turbulence. 

The next obvious question is, what kind of modifications 
does linear forcing need in order to reproduce aspects of 
a generic decaying turbulence, associated with (1591) and 
an evolution law L ~ i 7 ? Two immediate guesses are to 
consider a time-dependent rate A — A(t) or, to consider 
a time-dependent box whose size I evolves according to 
I ~ t 1 . The analysis of such possibilities is left for future 
work. 
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